1 Libraries and what not

## Base RF model
set.seed(42)

library(ggplot2)  # for autoplot() generic
library(dplyr)
library(sf)
library(stars)
library(tmap)

2 Data preprocessing

2.1 Glaciers

dat_sf <- st_read("./data/glacier_clim.shp")
## Reading layer `glacier_clim' from data source `/Users/u0784726/Dropbox/Data/devtools/glacier_mb/data/glacier_clim.shp' using driver `ESRI Shapefile'
## Simple feature collection with 95086 features and 31 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 67.47845 ymin: 27.4918 xmax: 103.8915 ymax: 45.35089
## Geodetic CRS:  GCS_unknown

2.2 Basemap

Load and crop basemap

basemap <- read_stars("~/Dropbox/Data/summer/natural_earth/HYP_50M_SR_W/HYP_50M_SR_W.tif")

basemap2 <- st_crop(basemap, st_bbox(dat_sf))
basemap2 <- st_rgb(basemap2)

3 ERA Data

3.1 t2m

t2m <- read_stars("~/Dropbox/Data/summer/era5/t2m_2001_2020_amjjas.nc")
st_crs(t2m) <- 4326
t2m <- st_crop(t2m, st_bbox(dat_sf))
t2m <- st_as_stars(t2m)

Make average 2006-2010

t2m_2008 <- t2m %>% 
  slice("time", 6:10) %>%
  st_apply(MARGIN = 1:2, FUN = mean)

Make average 2016-2020

t2m_2018 <- t2m %>% 
  slice("time", 16:20) %>%
  st_apply(MARGIN = 1:2, FUN = mean)

Make delta

t2m_d <- t2m_2018 - t2m_2008

3.2 tp

tp <- read_stars("~/Dropbox/Data/summer/era5/tp_2001_2020_ann.nc")
st_crs(tp) <- 4326
tp <- st_crop(tp, st_bbox(dat_sf))
tp <- st_as_stars(tp)

Make average 2006-2010

tp_2008 <- tp %>% 
  slice("time", 6:10) %>%
  st_apply(MARGIN = 1:2, FUN = mean)

Make average 2016-2020

tp_2018 <- tp %>% 
  slice("time", 16:20) %>%
  st_apply(MARGIN = 1:2, FUN = mean)

Make delta

tp_d <- tp_2018 - tp_2008

4 Maps

4.1 t2m

x <- c(t2m_2008, t2m_2018)
names(x) <- c("t2m_2008", "t2m_2018")
m1 <- tm_shape(basemap2) +
  tm_raster(palette = "Greys", alpha = 0.5) +
  tm_shape(x) +
  tm_raster(palette = "-inferno",
            alpha = 0.75, 
            style = "quantile") +
  tm_layout(legend.outside = TRUE)
print(m1)

Plot \(\delta\)

m1 <- tm_shape(basemap2) +
  tm_raster(palette = "Greys", alpha = 0.5) +
  tm_shape(t2m_d) +
  tm_raster(alpha = 0.75, 
            style = "quantile") +
  tm_layout(legend.outside = TRUE)
print(m1)

Plot \(\delta\) with mass balance

m1 <- tm_shape(t2m_d) +
  tm_raster(alpha = 0.75, 
            style = "quantile") +
  tm_shape(dat_sf) + 
    tm_symbols(col = "mb_mwea", 
             size = 0.25, 
             alpha = 0.75, 
             style = "quantile") +
  tm_layout(legend.outside = TRUE)
print(m1)

4.2 tp

x <- c(tp_2008, tp_2018)
names(x) <- c("tp_2008", "tp_2018")
m1 <- tm_shape(basemap2) +
  tm_raster(palette = "Greys", alpha = 0.5) +
  tm_shape(x) +
  tm_raster(palette = "YlGnBu",
            alpha = 0.75, 
            style = "quantile") +
  tm_layout(legend.outside = TRUE)
print(m1)

Plot \(\delta\)

m1 <- tm_shape(basemap2) +
  tm_raster(palette = "Greys", alpha = 0.5) +
  tm_shape(tp_d) +
  tm_raster(alpha = 0.75, 
            style = "quantile") +
  tm_layout(legend.outside = TRUE)
print(m1)

Plot \(\delta\) with mass balance

m1 <- tm_shape(tp_d) +
  tm_raster(alpha = 0.75, 
            style = "quantile") +
  tm_shape(dat_sf) + 
    tm_symbols(col = "mb_mwea", 
             size = 0.25, 
             alpha = 0.75, 
             style = "quantile") +
  tm_layout(legend.outside = TRUE)
print(m1)